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Abstract. A popular theory of self-organized criticality relates the critical behavior of 
driven dissipative systems to that of systems with conservation. In particular, this theory 
predicts that the stationary density of the abelian sandpile model should be equal to the 
threshold density of the corresponding fixed-energy sandpile. This "density conjecture" has 
been proved for the underlying graph Z. We show (by simulation or by proof) that the 
density conjecture is false when the underlying graph is any of 1?, the complete graph Km 
the Cayley tree, the ladder graph, the bracelet graph, or the flower graph. Driven dissipative 
sandpiles continue to evolve even after a constant fraction of the sand has been lost at the 
sink. These results cast doubt on the validity of using fixed-energy sandpiles to explore the 
critical behavior of the abelian sandpile model at stationarity. 



1. Introduction 

In a widely cited series of papers |DVZ98t IVDMZ981 IDMVZnOi IVOMZOOl IMPFS+Olj . 
Dickman, Muiioz, Vespignani and Zapperi (DMVZ) developed a theory of self-organized 
criticality as a relationship between driven dissipative systems and systems with conservation. 
This theory predicts a specific relationship between the classical abelian sandpile model of 
Bak, Tang, and Wiesenfeld |BTW87] . a driven system in which particles added at random 
dissipate across the boundary, and the corresponding "fixed-energy sandpile," a closed system 
in which the total number of particles is conserved. 

In this introduction, we briefly define these two models and explain the conjectured rela- 
tionship between them in the DMVZ paradigm of self-organized criticality. In particular, we 
focus on the prediction that the stationary density of the driven dissipative model equals the 
threshold density of the fixed-energy sandpile model. In lsection 21 we present data from large- 
scale simulations which strongly indicate that this conjecture is false on the two-dimensional 
square lattice Z^. In the subsequent sections we expand on the results announced in [FLWIO] 
by examining the conjecture on some simpler families of graphs in which we can provably 
refute it. 

The difference between the stationary and threshold densities on most of these graphs is 
fairly small — typically on the order of 0.01% to 0.2% — which explains why many previous 
simulations did not uncover them. The exception is the early experiments by Grassberger 
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and Manna |GM90j . who clearly identified this discrepancy, at least in dimensions 4 and 
higher. Later studies focused on dimension 2 and missed the discrepancy. 

In some more recent papers such as |BM08] . the DMVZ paradigm is explicitly restricted to 
stochastic models. In other recent papers |VD05t IdCVdSDOQ] it is claimed to apply both to 
stochastic and deterministic sandpiles, although these papers focus on stochastic sandpiles, 
for the reason that deterministic sandpiles are said to belong to a different universality class. 

Despite our contrary findings, we believe that the central idea of the DMVZ paradigm is a 
good one: even in the deterministic case, the dynamics of a driven dissipative system should 
in some way reflect the dynamics of the corresponding conservative system. Our results point 
to a somewhat different relationship than that posited in the DMVZ series of papers: the 
driven dissipative model exhibits a second-order phase transition at the threshold density of 
the conservative model. We explain this transition in [section 3l 

Bak, Tang, and Wiesenfeld jBTW87] introduced the abelian sandpile as a model of self- 
organized criticality; for mathematical background, see |Red05] . The model begins with a 
collection of particles on the vertices of a finite graph. A vertex having at least as many 
particles as its degree topples by sending one particle along each incident edge. A subset of 
the vertices are distinguished as sinks: they absorb particles but never topple. A single time 
step consists of adding one particle at a random site, and then performing topplings until 
each non-sink vertex has fewer particles than its degree. The order of topplings does not 
affect the outcome [DhaQOj . The set of topplings that occur before the system stabilizes is 
called an avalanche. 

Avalanches can be decomposed into a sequence of "waves" so that each site topples at 
most once during each wave. Over time, sandpiles evolve toward a stationary state in which 
the waves exhibit power-law statistics |KLGPOO] (though the full avalanches seem to exhibit 
multifractal behavior |DMST98| IKMS05] ). Power law statistics are a hallmark of criticality, 
and since the stationary state is reached apparently without tuning of a parameter, the 
model is said to be self- organized critical. 

To explain how the sandpile model self-organizes to reach the critical state, Dickman et 
al. |DVZ98t IDMVZOO] introduced an argument which soon became widely accepted: see, 
for example, |Sor041 Ch. 15.4.5] and |MQ05i iFdBROSl [RS09] . Despite the apparent lack of 
a free parameter, they argued, the dynamics implicitly involve the tuning of a parameter 
to a value where a phase transition takes place. The phase transition is between an active 
state, where topplings take place, and a quiescent "absorbing" state where topplings have 
died out. The parameter is the density, the average number of particles per site. When 
the system is quiescent, addition of new particles increases the density. When the system 
is active, particles are lost to the sinks via toppling, decreasing the density. The dynamical 
rule "add a particle when all activity has died out" ensures that these two density changing 
mechanisms balance one another out, driving the system to the threshold of instability. 

To explore this idea, DMVZ introduced the fixed- energy sandpile model (FES), which 
involves an explicit free parameter (, the density of particles. On a graph with vertices, the 
system starts with (N particles at vertices chosen independently and uniformly at random. 
Unlike the driven dissipative sandpile described above, there are no sinks and no addition 
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of particles, so the total number of particles is conserved. Subsequently the system evolves 
through toppling of unstable sites. Usually the parallel toppling order is chosen: at each time 
step, all unstable sites topple simultaneously. In the mathematical literature, this system 
goes by the name of parallel chip-firing |BG92t IBLS91j . Toppling may persist forever, or it 
may stop after some finite time. In the latter case, we say that the system stabilizes; in the 
terminology of DMVZ, it reaches an "absorbing state." 

A common choice of underlying graph is the n x n square grid with periodic boundary 
conditions. It is believed, and supported by simulations |BCFV03] . that there is a threshold 
density such that for ( < the system stabilizes with probability tending to 1 as n — )■ oo; 
and for ( > (c, with probability tending to 1 the system does not stabilize. 

For the driven dissipative sandpile on the n x n square grid, as n — t- oo the stationary 
measure has an infinite- volume limit [A J04j . which is a measure on sandpiles on 7?. It 
turns out that one gets the same limiting measure whether the grid has periodic or open 
boundary conditions, and whether there is one sink vertex or the whole boundary serves 
as a sink |AJ04] (see also [PemQlj for the corresponding result on random spanning trees). 
The statistical properties of this limiting measure have been much studied |Pri94a[ IJPROGj . 
Of particular interest is the stationary density (g of Z^, defined as the expected number of 
particles at a fixed site. Grassberger conjectured that (s is exactly 17/8, and it is now known 
that C. = 17/8 ± 10-12 [JPK06]. 

DMVZ believed that the combination of driving and dissipation in the classical abelian 
sandpile model should push it toward the critical density (c of the fixed-energy sandpile. 
This leads to a specific testable prediction, which we call the Density Conjecture. 

Conjecture 1.1 (Density Conjecture, |VDMZOO] ). On the square grid, = 17/8. 

One can also formulate a density conjecture for more general graphs, where it takes the 
form Cc = (s- We give precise definitions of the densities (c and (s in section [2l 

Vespignani et al. [VDMZOO] write of the fixed-energy sandpile on the square grid, "the 
system turns out to be critical only for a particular value of the energy density equal to that 
of the stationary, slowly driven sandpile." They add that the threshold density Cc of the 
fixed-energy sandpile is "the only possible stationary value for the energy density" of the 
driven dissipative model. In simulations they find Cc = 2.1250(5), adding in a footnote "It 
is likely that, in fact, 17/8 is the exact result." 

Munoz et al. [MDPS"'"Ol] have also expressed this view, asserting that "FES are found to 
be critical only for a particular value C = Cc (which as we will show turns out to be identical 
to the stationary energy density of its driven dissipative counterpart)." 

Our goal in the present paper is to demonstrate that the density conjecture is more prob- 
lematic than it first appears. Table [T] presents data from large-scale simulations which 
strongly suggest that Cc is close to but not exactly equal to 17/8 (see also Table [3]). We 
further consider several other families of graphs, including some for which we can determine 
the exact values Cc and Cs analytically. We find that they are close, but not equal. 

All now known information on the threshold density Cc and stationary density Cs is sum- 
marized in Table [2l The only graph on which the two densities are known to be equal is Z 
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n 



64 1 28 2 5 6 5 1 2 1024 2048 4096 8192 16384 



Table 1. Summary of our fixed-energy sandpile simulations on n x n tori Z^, 
giving our empirical estimate of the threshold density Cc(^^)- The standard 
deviation in each of our estimates of Cc(^^) is 4 x 10^''. The data are well 
approximated by Cc(ZD = 2.1252881 ± 3 x 10"^ - (0.390 ± 0.001)n-i-^ as 
shown in the graph. (The error bars are too small to be visible, so the data 
are shown as points.) We conclude that the asymptotic threshold density 
Cc(Z^) is 2.125288 to six decimal places. In contrast, the stationary density 
Cil?) is 2.125000000000 to twelve decimal places. 
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Table 2. Stationary and threshold densities for different graphs. Exact values 
are in bold. 



|MQ05 , rFdBR05[ IFdBMR09] . On all other graphs we examined, with the possible exception 
of the 3-regular tree, it appears that Cc 7^ C- 

Taken together, these results show that the conclusions of [MDPS"'"Ol] that "FES are 
shown to exhibit an absorbing state transition with critical properties coinciding with those 
of the corresponding sandpile model" deserve to be re-evaluated. One hope of the DMVZ 
paradigm was that critical features of the driven dissipative model, such as the exponents 
governing the distribution of avalanche sizes and decay of correlations, might be more easily 
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studied in FES by examining the scaling behavior of these observables as C t Cc- However, 
the failure of the density conjecture suggests that the two models may not share the same 
critical features. 

As Grassberger and Manna observed |GM90j . the value of the FES threshold density 
depends on the choice of initial condition. One might consider a more general version of FES, 
namely adding (C — Co)-^ particles at random to a "background" configuration r of density Co 
already present on the grid. For example, taking r to be the deterministic configuration on 
oild — l particles everywhere, by |FLP10l Prop 1.4] we obtain a threshold density of Cc = 
Id — I. Many interesting questions present themselves: for instance, for which background 
does Cc take the smallest value, and for which backgrounds do we obtain Cc = Cs? It would 
also be interesting to replicate the phase transition for driven sandpiles (see [section 3"!) for 
different background configurations. 



2. Sandpiles on the square grid 1? 

In this section we give precise definitions of the stationary and threshold densities, and 
present the results of large-scale simulations on J? . The definitions in this section apply 
to general graphs, but we defer the discussion of results about other graphs to subsequent 
sections. 

2.1. The driven dissipative sandpile and the stationary density Cs- Let G = (V, E) 

be a finite graph, which may have loops and multiple edges. Let S (Z V he a nonempty set of 
vertices, which we will call sinks. The presence of sinks distinguishes the driven dissipative 
sandpile from its fixed-energy counterpart. To highlight this distinction, throughout the 
paper, graphs denoted with a "hat" as in G have sinks, and those without a hat as in G do 
not. 

For vertices v,w E V, write a^,^ = aw,v for the number of edges connecting v and w, and 



for the number of edges incident to v. A sandpile (or "configuration") rj on G is a map 

r]:V Z>o. 

We interpret t]{v) as the number of sand particles at the vertex v; we will sometimes call 
this number the height of v in t]. 

A vertex v ^ S is called unstable if 77(f) > d^. An unstable vertex can topple by sending 
one particle along each edge incident to v. Thus, toppling v results in a new sandpile rj' 
given by 

7]' = r] + 

where 

v,w, V 
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Figure 1. The square grid I?. 



Sinks by definition are always stable, and never topple. If all vertices are stable, we say 
that 7] is stable. 

Note that toppling a vertex may cause some of its neighbors to become unstable. The 
stabilization r]° of is a sandpile resulting from toppling unstable vertices in sequence, until 
all vertices are stable. By the abelian property |Dha90] . the stabilization is unique: it does 
not depend on the toppling sequence. Moreover, the number of times a given vertex topples 
does not depend on the toppling sequence. 

The most commonly studied example is the n x n square grid graph, with the boundary 
sites serving as sinks ([Figure ij). 



The driven dissipative sandpile model is a continuous time Markov chain {rit)t>o whose 
state space is the set of stable sandpiles on G. Let V = V \ S he the set of vertices that 
are not sinks. At each site v E V, particles are added at rate 1. When a particle is added, 
topplings occur instantaneously to stabilize the sandpile. Writing (Jt{v) for the total number 
of particles added at v before time t, we have by the abelian property 

Vt = {cTtT. 

Note that for fixed t, the random variables (Tt{v) for v E V are independent and have the 
Poisson distribution with mean t. 

The model just described is most commonly known as the abelian sandpile model (ASM), 
but we prefer the term "driven dissipative" to distinguish it from the fixed-energy sandpile 
described below, which is also a form of ASM. "Driven" refers to the addition of particles, 
and "dissipative" to the loss of particles absorbed by the sinks. 

Dhar |Dha90] developed the burning algorithm to characterize the recurrent sandpile 
states, that is, those sandpiles rj for which, regardless of the initial state, 

Pr(?7j = rj for some t) = 1. 

Lemma 2.1 (Burning Algorithm |Dha90j ). A sandpile rj is recurrent if and only if every 
non-sink vertex topples exactly once during the stabilization of rj + Ylses^s, where the sum 
is over sink vertices S . 



THE APPROACH TO CRITICALITY IN SANDPILES 



7 



The recurrent states form an abelian group under the operation of addition followed by 
stabilization. In particular, the stationary distribution of the Markov chain 77^ is uniform on 
the set of recurrent states. 

The combination of driving and dissipation organizes the system into a critical state. To 
measure the density of particles in this state, we define the stationary density Cs{G) as 



Cs{G) = E, 



where V = V\S , and fi is the uniform measure on recurrent sandpiles on G. The stationary 
density has another expression in terms of the Tutte polynomial of the graph obtained from G 
by collapsing the set S of sinks to a single vertex; see [section 41 

Most of the graphs we will study arise naturally as finite subsets of a locally finite graph 
r, i.e., r is a countably infinite graph in which every vertex has finite degree. (We also study 
the complete graph and the flower graph, which do not arise in this way.) Let Gn for ri > 1 
be a nested family of finite induced subgraphs with [J G„ = F. As sinks in G„ we take the 
set of boundary vertices 

Sn — Gn Gn-l- 

In cases where the free and wired limits are different, such as on regular trees, we will choose 
a sequence Gn corresponding to the wired limit. We denote by Hn the uniform measure on 
recurrent configurations on G„. 

We are interested in the stationary density 

C(r) := lim C{Gn). 

n— >oo 

When r = Z*^, it is known that the infinite-volume limit of measures /i = lim^^oo /^n exists 
and is translation-invariant |AJ04] . In this case it follows that the limit defining Cs(r) exists 
and equals 

Cs = E,[r]{0)]. 

For other families of graphs we consider, we will show that the limit defining Cs(r) exists. 
Much is known about the limiting measure fi in the case F = Z^. The following expressions 

have been obtained for (s and the single site height probabilities. The symbol = denotes 
expressions that are rigorous up to a conjecture |JPR06j that a certain integral, numerically 
evaluated as 0.5 ± 10~^^, is exactly 1/2. 

C(Z^) = 17/8 |JPR06j . 

^{r^(x) = o} = ^-^ mm, 

/xMx) = l} = i-^-^ + l| |PH94^[JPR06], 
/i{r/(x) = 2} = | + i-i| [JPR06] , and 
^{r^(x) = 3} = l-^ + ^ + ^ [JPR06]. 
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2.2. The fixed-energy sandpile model and the threshold density (c- Next we describe 
the fixed-energy sandpile model, in which the driving and dissipation are absent, and the 
total number of particles is conserved. As before, let G be a finite graph, possibly with loops 
and multiple edges. Unlike the driven dissipative model, we do not single out any vertices 
as sinks. The fixed-energy sandpile evolves in discrete time: at each time step, all unstable 
vertices topple once in parallel. Thus the configuration rjj+i at time j -|- 1 is given by 

Vj+i = Vj+Yl 

where 

Uj = {v e V : rjjiy) > d^} 

is the set of vertices that are unstable at time j. We say that t]o stabilizes if toppling 
eventually stops, i.e. Uj = for all sufficiently large j. 

If rjo stabilizes, then there is some site that never topples |Tar88j (see also |FdBMR09| 
Theorem 2.8, item 4] and [FLPlOl Lemma 2.2]). Otherwise, for each site x, let j{x) be 
the last time x topples. Choose a site x minimizing j{x). Then each neighbor y of x has 
j{y) > so y topples at least once at or after time j{x). Thus x receives at least 

additional particles and must topple again after time j{x), a contradiction. This criterion is 
very useful in simulations: as soon as every site has toppled at least once, we know that the 
system will not stabilize. 

Let {crx{v))x>o be a collection of independent Poisson point processes of intensity 1, indexed 
by the vertices of G. So each (Tx{v) has the Poisson distribution with mean A. We define the 
threshold density of G as 

UG) = EAe, 

where 

Ac = sup{A : (Xx stabilizes}. 

We expect that Ac is tightly concentrated around its mean when G is large. Indeed, if F is an 
infinite vertex-transitive graph, then the event that ax stabilizes on F is translation- invariant. 
By the ergodicity of the Poisson product measure, this event has probability or 1. Since 
this probability is monotone in A, there is a (deterministic) threshold density Cc(r), such 
that 

Prlax stabilizes on Fl = < ' "^"^^ ^ 

^ J \o, A>Cc(r). 

We expect the threshold densities on natural families of finite graphs to satisfy a law of large 
numbers such as the following. 

Conjecture 2.2. With probability 1, 

Ac(Z^) -> Cc(Z^) as n oo. 

Previous simulations (n = 160 |DVZ98j : n = 1280 |VDMZ98] ) to estimate the threshold 

density Cc(Z^) found a value of 2.125, in agreement with the stationary density Cs(Z^) = 17/8. 
By performing larger-scale simulations, however, we find that (c exceeds (s- 
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grid size 


^samples 


UK) 


distribution of height h of sand 


#topphngs 




Pr[/i = 0] 


Pr[/i = 1] 


Pr[/i = 2] 


Pr[/i = 3] 






268435456 


2.124956 


0.073555 


0.173966 


0.306447 


0.446032 


0.197110 


128^ 


67108864 


2.125185 


0.073505 


0.173866 


0.306567 


0.446062 


0.197808 


256^ 


16777216 


2.125257 


0.073488 


0.173835 


0.306609 


0.446068 


0.198789 


512^ 


4194304 


2.125279 


0.073481 


0.173826 


0.306626 


0.446067 


0.200162 


10242 


1048576 


2.125285 


0.073479 


0.173822 


0.306633 


0.446066 


0.201745 


20482 


262144 


2.125288 


0.073478 


0.173821 


0.306635 


0.446065 


0.203378 


40962 


65536 


2.125288 


0.073477 


0.173821 


0.306637 


0.446064 


0.205323 


81922 


16384 


2.125288 


0.073477 


0.173821 


0.306638 


0.446064 


0.206475 


163842 


4096 


2.125288 


0.073478 


0.173821 


0.306638 


0.446064 


0.208079 


Z'^ (stationary) 


2.125000 


0.073636 


0.173900 


0.306291 


0.446172 





Table 3. Fixed-energy sandpile simulations on n x n tori Z2. The third 
column gives our empirical estimate of the threshold density W^n)- The next 
four columns give the empirical distribution of the height of a fixed vertex in 
the stabilization {crx)°, for A just below Ac. Each estimate of the expectation 
CcC^n) of marginals Pr[/i = i] has standard deviation less than 4- 10"''. 
The total number of topplings needed to stabilize a\ appears to scale as n^. 



Table [3] summarizes the results of our simulations, which indicate that Cc(^^) equals 
2.125288 to six decimal places. In each random trial, we add particles one at a time at 
uniformly random sites of the n x n torus. After each addition, we perform topplings until 
either all sites are stable, or every site has toppled at least once since the last addition. For 
deterministic sandpiles on a connected graph, if every site topples at least once, the system 
will never stabilize |FdBMR09| IFLPlOj . We record m/v? as an empirical estimate of the 
threshold density, where m is the maximum number of particles for which the configuration 
stabilizes. We then average these empirical estimates were over many independent trials. 
The one-site marginals we report are obtained from the stable configuration just before the 
(m + 1)^* particle was added, and the number of topplings reported is the total number of 
topplings required to stabilize the first m particles. 

We used a random number generator based on the Advanced Encryption Standard (AES- 
256), which has been found to exhibit good randomness properties. Our simulations were 
conducted on a High Performance Computing (HPC) cluster of computers. 

3. Sandpiles on the bracelet 

Next we examine a family of graphs for which we can determine and exactly and 
prove that they are not equal. Despite this inequality, we show that an interesting connection 
remains between the driven dissipative and conservative dynamics: the threshold density of 
the conservative model is the point at which the driven dissipative model begins to lose a 
macroscopic amount of sand to the sink. 



The bracelet graph Bn (Figure 2) is a multigraph with vertex set Z„ (the n-cycle) with the 



usual edge set + 1 mod n) : Q <i <n} doubled. Thus all vertices have degree 4. The 
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graph Bn is the same, except that vertex is distinguished as a sink from which particles 
disappear from the system. We denote by B^o the infinite path Z with doubled edges. 

For A > 0, let a\ be the configuration with Poisson(A) particles independently on each 
site of Bn- Let rjx = {crx)° be the stabilization of ax, and let 

^ n—l 
a:=l 

be the final density. The following theorem gives the threshold and stationary densities of 
the infinite bracelet graph B^o, and identifies the ri — )■ oo limit of the final density p„(A) as 
a function of the initial density A. 

Theorem 3.1. For the bracelet graph, 

(1) The threshold density Cc(-Boo) is the unique positive root of( = | — |e~^'' (numerically, 
Cc = 2.496608;. 

(2) The stationary density Cs(-Boo) is 5/2. 

(3) p„(A) — )■ p(A) in probability as n ^ oo, where 



p(A) = min y\, 

Part 3 of this theorem shows that the final density undergoes a second-order phase tran- 




sition at (c- the derivative of p(A) is discontinuous at A = Cc ( [Figure 3"| ). Thus in spite of 
the fact that (s 7^ Cc, there remains a connection between the conservative dynamics used 
to define Cc and the driven-dissipative dynamics used to define Cs- For A < Cc, very little 
dissipation takes place, so the final density equals the initial density A; for A > Cc a sub- 
stantial amount of dissipation takes place, many particles are lost to the sink, and the final 
density is strictly less than the initial density. The sandpile continues to evolve as A increases 
beyond Cc! in particular its density keeps changing. 

We believe that this phenomenon is widespread. As evidence, in [section 5) we introduce 
the "flower graph," which looks very different from the bracelet, and prove (in [Theorem 5.5|) 
that a similar phase transition takes place there. 

For the proof of [Theorem 3.11 we compare the dynamics of pairs of particles on the bracelet 
graph to single particles on Z. At each vertex x of the bracelet, we group the particles starting 




Figure 2. The bracelet graph i?2o- 
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1 2 Cc 3 2.4 Cc2.5 2.6 

Figure 3. Density p(A) of the final stable configuration as a function of initial 
density A, for the driven sandpile on the bracelet graph Bn as n — )■ cxo. A 
second-order phase transition occurs at A = (c- Beyond this transition, the 
density of the driven sandpile continues to increase, approaching the stationary 
density Cs from below. 



at X into pairs, with one "passive" particle left over if ax{x) is odd. Since all edges in the 
bracelet are doubled, we can ensure that in each toppling the two particles comprising a pair 
always move to the same neighbor, and that the passive particles never move. The toppling 
dynamics of the pairs are equivalent to the usual abelian sandpile dynamics on Z. 
We recall the relevant facts about one- dimensional sandpile dynamics: 

• In any recurrent configuration on a finite interval of Z, every site has height 1, except 
for at most one site of height 0. Therefore, C = 1 |Red05j . 

• On Z, an initial configuration distributed according to a nontrivial product measure 
with mean A stabilizes almost surely (every site topples only finitely many times) if 
A < 1, while it almost surely does not stabilize (every site topples infinitely often) if 
A > 1 |FdBMR09j . Thus, Cc = 1- 

Proof of I Theorem 3. 1\ parts 1 and 2. Given A > 0, let A* be the pair density ELcrA(a;)/2j , 
and let 

\2m+l 1 

i.odd(A) = e-^5^-^-- = ^(l-e-^^). 
^-^ (Zm +1)1 2 

be the probability that a Poisson(A) random variable is odd. Then A and A* are related by 

A = 2A*+Podd(A). (1) 

The configuration a\ stabilizes on Boo if and only if the pair configuration (y\ stabilizes on Z. 
Thus Cc(5oo)* = Cc(Z)- Setting A = Cc{Boo) in dlD, using the fact that Cc(Z) = 1, and that A* 
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is an increasing function of A > 0, we conclude that Cc(-Boo) is the unique positive root of 

C = 2+Podd(C), 

or C = I — |e~^^. This proves part 1. 

For part 2, by the burning algorithm, a configuration a on is recurrent if and only 
if it has at most one site with fewer than two particles. Thus, in the uniform measure on 
recurrent configurations on i?^, 

Pr(cr(x) = 2) = Pr(a(x) = 3) = ^ - ^, Pr(a(x) = 0) = Vy{o{x) = 1) = ^. 

We conclude that Cs(-Sn) = IEo-(a;) = | - | | as n oo. □ 

To prove part 3 of [Theorem 3.1t we use the following lemma, whose proof is deferred to 
the end of this section. Let be the n-cycle with vertex distinguished as a sink. Let 
be a sandpile on Z„ distributed according to a product measure (not necessarily Poisson) 
of mean A. Let r/^ be the stabilization of cr^, and let p^(A) = ■^^^YTx=\'^\^^) ^^^^ 
density after stabilization. 

Lemma 3.2. On Z„, we have p'„(A) — i- min(A, 1) in 'probability. 

Proof of \Theorem 3.11 part 3. Let rjx be the stabilization of ax on Bn, and let ri\ be the 
stabilization of = [_o'x/2\ on Then 

r]x{x) = 2r]l{x) + UJx{x) (2) 

where ux{x) = (Tx{x) — 2(tI{x) is 1 or accordingly as o"a(x) is odd or even. Let 

n-1 



x=l 

be the final density after stabilization of cr^ on Z„. Then 

n—l 

p„(A) = 2p;(A) + — -Y.^^(' 

n—l ^ 



x=l 

By the weak law of large numbers, Y12=i ^>^i^) ~^ Podd(A) in probability as n — )■ oo. If 
A < Co then A* < 1, so by lLemma 3.2[ p*(A) — > A* in probability, and hence 

p„(A) ^2A*+Podd(A) = A 

in probability. If A > then A* > 1, so by lLemma 3.21 p*(A) — 1 in probability, hence 

5 - e"^-^ 

p„(A) 2 +f»odd(A) = 

in probability. This proves part 3. □ 
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Proof of \Lemma 3.21 We may view Z„ as the path in Z from a„ = — [ra/2j to = \n/2], 
with both endpoints serving as sinks. For x G Z„, let Un{x) be the number of times that 
X topples during stabilization of the configuration cr^ on Z„. Let Uoo{x) be the number 
of times x topples during stabilization of cr^ on Z. The procedure of "toppling in nested 
volumes" |FdBMR09] shows that Un{x) t Uoo{x) as n — >■ oo. 

We consider first A < 1. In this case Uoo{x) is finite almost surely (a.s.). The total number 
of particles lost to the sinks on Z„ is 'U„(a„ + 1) + m„(6„ — 1), so the final density is given by 

bn-l 

O-xix) - Un{an + 1) " Un{bn " 1) 



n — 1 



-1 



By the law of large numbers, ^ o-x^x) — )■ A in probability as n — > oo. Since Uoo{x) is a.s. 
finite, we have +^^+""(^"~^) _i. g probability, so p^(A) — )■ A in probability. 

Next we consider A > 1. In this case we have t Uoo{x) = oo, a.s. Let p{n,x) = 

Pr('u„(x) = 0) be the probability that x G Z„ does not topple. By the abelian property, 
adding sinks can not increase the number of topplings, so 

p{n, x) < p{m, 0) 



where m = min(x — an,bn — x). Let 



x=an+l 

be the number of sites in Z^ that do not topple. Since UniO) t c>o a.s., we have p(n, 0) I 0, 
hence 

E— = - V p(n,a;) < - Vp(m,0) ^ 
n n ^-^ n ^-^ 

i'=an+l m=l 

as n — )■ oo. Since l^n > it follows that Ynjn — ?■ in probability. 

In an interval where every site toppled, there can be at most one empty site. We have 
y„ + 1 such intervals. Therefore, the number of empty sites is at most 21^ + 1. Hence 

n - 2r„ - 2 , , , , 

< P;(A) < 1. 

n — 1 

The left side tends to 1 in probability, which completes the proof. □ 

4. Sandpiles on the complete graph 

Let Kn be the complete graph on n vertices: every pair of distinct vertices is connected by 
an edge. In Kn., one vertex is distinguished as the sink. The maximal stable configuration 
on Kn has density n— 2, while the minimal recurrent configurations have exactly one vertex of 
each height 0, 1, . . . , n — 2, hence density The following result shows that the stationary 
and threshold densities are quite far apart: is close to the minimal recurrent density, 
while Cc is close to the maximal stable density. 
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Theorem 4.1. 

C(ir„) = - + o(V^) 

Cc{Kn) >n- Oi^/nAogn). 

The proof uses an expression for the stationary density (s in terms of the Tutte polynomial, 
due to C. Merino Lopez |ML97j . Our application will be to the complete graph, but we state 
Merino Lopez' theorem in full generality. Let G = {V, E) be a connected undirected graph 
with n vertices and m edges. Let v be any vertex of G, and write G for the graph G with v 
distinguished as a sink. Let d be the degree of v. 

Recall that the Tutte polynomial To{x,y) is defined by 



TG(x,t/) = 5^(x-ir(^)-^(^)(y 



ACE 

where c{A) denotes the number of connected components of the spanning subgraph {V,A). 
Theorem 4.2 ( |ML97] ). The Tutte polynomial TG{x,y) evaluated at x = 1 is given by 

Tg(1,i/) = /-"5^1/I'^I 

where the sum is over all recurrent sandpile configurations a on G, and \a\ denotes the 
number of particles in a. 

Differentiating and evaluating at ?/ = 1, we obtain 

Y.{d-m+\a\). (3) 



y=l 



Referring to the definition of the Tutte polynomial, we see that Tg(1, 1) is the number 
of spanning trees of G, and that the left side of ([3]) is the number of spanning unicyclic 
subgraphs of G. (In evaluating at x = ?/ = 1, we interpret 0° as 1.) The number of 
recurrent configurations equals the number of spanning trees of G, so the stationary density 
C,s niay be expressed as 



Combining these expressions yields the following: 
Corollary 4.3. 

n \ i^\G) ^ 

where k{G) is the number of spanning trees of G, and u{G) is the number of spanning 
unicyclic subgraphs ofG. 
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Note that m — d is the minimum number of particles in a recurrent configuration, so the 
ratio u{G) / k{G) can be interpreted as the average number of excess particles in a recurrent 
configuration. 

Everything so far applies to general connected graphs G. The following is specific for the 
complete graph. 

Theorem 4.4 (Wright |Wri77] ). The number of spanning unicyclic subgraphs of K„ 



IS 



+ o(l) 



Proof of \Theorem 4 -11 For Kn we have 

n(n — 1) 

m — a = 



(n-l) 



2 ^ ' 2 

From Corollary 4.3 [Theorem 4.4] and Cayley's formula hi{Kn 



n - 2){n - 1) 

we obtain 



n 



riK\ 1 ( {n-2){n-\) , u{K^) 
^'^""-^ = n [ 2 + ^) 



n 



+ 



On the other hand, if we let 



A = n — 2 y 72 log n 

and start with a{v) ~Poisson(A) particles at each vertex v of Kn, then for all v 



Pr[cr('u) >n]<—. 



So 



Pr[cr(f) > n for some v] < 

n 



in other words, with high probability no topplings occur at all. Thus 

Pr (^Ac{Kn) >n- 2^n lognj > 1 - - 
which completes the proof. 



□ 



One might guess that the large gap between Cs and Cc is related to the small diameter 
of Kn'- since the sink is adjacent to every vertex, its effect is felt with each and every toppling. 
This intuition is misleading, however, as shown by the lollipop graph L„ consisting of Kn 
connected to a path of length n, with the sink at the far end of the path. Since L„ has the 



same number of spanning trees and unicyclic subgraphs as Kn-, we have by Corollary 4.3 

I m \ 
I , " ^ ^ 

n{n - 1) , ^ , u[Ln) 



Cs{Ln) — — 

2n 



+ n-l 



KiLr. 



\ 



n 



4+0(v^). 
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On the other hand, by first stabihzing the vertices on the path, close to half of which end 
up in the sink without reaching the Kn-, it is easy to see that with high probability 

2n 

Ac(^n) > — -0{^Jn\ogn). 

5. SaNDPILES on THE FLOWER GRAPH 

An interesting feature of parallel chip-firing is that further phase transitions appear above 
the threshold density Cc- On a finite graph G = {V,E), since the time evolution is deter- 
ministic, the system will eventually reach a periodic orbit: for some positive integer m, we 
have r]t+m = Vt for all sufficiently large t. The activity density, pa, measures the proportion 
of vertices that topple in an average time step: 

1 1 

Pa{X) = Ea lim -J2^Y1 Mvsi^)>d.}- 
s=o ^ xev 

The expectation Ea refers to the initial state rjo, which we take to be distributed according 
to the Poisson product measure with mean A. Note that the limit in the definition of pa can 
also be expressed as a finite average, due to the eventual periodicity of the dynamics. 

Bagnoli et al. |BCFV03] observed that pa tends to increase with A in a sequence of flat 
steps punctuated by sudden jumps. This "devil's staircase" phenomenon is so far explained 
only on the complete graph [LevOSj : The number of flat stairs increases with n, and in the 
n — >■ (X) limit there is a stair at each rational number height pa = p/q- 

On the cycle Z„ |Dal06j there are just two jumps: at A = 1, the activity density jumps from 
to 1/2, and at A = 2, from 1/2 to 1. For the n x n torus, simulations |BCFV03] indicate 
a devil's staircase, which is still not completely understood despite much effort [CDVVOG] . 

In this section we study the "flower" graph, which was designed with parallel chip-firing 
in mind: the idea is that a graph with only short cycles should give rise to short period 
orbits under the parallel chip-firing dynamics. We find that there are four activity density 
jumps fITheorem 5.41) . In addition, we determine the stationary and threshold densities of 
the flower graph, and find a second-order phase transition at (c fITheorem 5.5|) . 



The flower graph consists of a central site together with n > 1 petals (Figure 4). Each 



petal consists of two sites connected by an edge, each connected to the central site by an 
edge. Thus the central site has degree 2n, and all other sites have degree 2. The number of 
sites is 2n + 1. The graph Fn is the same, except one petal serves as sink. 

Recall that we defined the density of a configuration as the total number of particles, 
divided by the total number of sites. Since the flower graph is not regular, the central site 
has a different expected number of particles than the petal sites. 

Proposition 5.1. For parallel chip-firing on the flower graph Fn, every configuration has 
eventual period at most 3. 

The proof uses the following two lemmas. 

Lemma 5.2 ( |Pri94bj Lemma 2.5). // the eventual period is not 1, then after some finite 
time, every site x has height at most 2dx — 1. 
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We also use an observation from |BCFV03] : it was stated and proved there for Z^, but the 
same proof works for general graphs. 

Lemma 5.3 ( |BCFVd3] ). Let two height configurations rj and ^ be "mirror images" of each 
other, that is, r]{x) = 2dx — 1 — ^(x) for all x. Then after performing for each a parallel 
chip-firing time step, the two configurations are again mirror images of each other. 



Proof of Proposition 5.1 Suppose at time t the model has settled into periodic orbit, and 
the period is not 1. Then by ILemma 5.2\ at time t every petal site has height at most 3. 
Say that a petal is in state ij if it has i particles at one site and j particles at the other site. 
A priori there are 10 possible petal states, listed below, where each one has two possible 
successor states, depending on whether or not the central site is stable. If a petal is in state 
ij, then by S{ij) we denote the state that it is in after one time step in which the central site 
does not topple, and likewise by U(zj) after one time step in which the central site topples. 

state S (state) U (state) 



00 


00 


11 


01 


01 


12 


02 


01 


12 


03 


11 


22 


11 


11 


22 


12 


02 


13 


13 


12 


23 


22 


11 


22 


23 


12 


23 


33 


22 


33 



From this we see that a petal will be in state 00 only if the central site is always stable, 
and consequently each site is always stable, in which case the period is 1. Similarly, petal 
state 33 only occurs if the central site is unstable each step, in which case each site must 
be unstable each step, and the period is again 1. State 03 is not a successor of any state of 
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these states, so it will not be a periodic petal state either. Thus the set of allowed periodic 
petal states is {01, 02, 11, 12, 13, 22, 23}. 

If the central site is stable every other time step, then the possible petal states are 12 — >■ 
02 ^ 12, 22 ^ 11 ^ 22, and 13 ^ 12 ^ 13, each of which has period 2. Then the period 
of the entire configuration is 2. 

Thus if the period is larger than 2, the central site must be stable for at least two con- 
secutive time steps, or else unstable for at least two consecutive time steps. We will label a 
time step S if the central site is stable in that time step, otherwise we label it U. So, if the 
period is larger than 2 we will see SS or UU in the time evolution. In the latter case, we can 
study the mirror image, which will have the same period, and for which we will see SS. 

Eventually the central site must be unstable again, since otherwise the period would be 1. 
Therefore, we can examine three time steps labeled SSU. Examining the evolution of the 
central site together with the petals, we see 



S 


S 


U 




01, 02 


01 


01 


12 


12, 23 


02 


01 


12 


11, 22 


11 


11 


22 


23 


12 


02 


12 



Whenever we have SSU, during the second and third time steps each petal contributes at 
most two particles to the central site, while the central site topples, so the central site must 
again be stable. Thus SSUU cannot occur, and we see SSUS. 

There are two cases for what the central site does next. Let us first consider SSUSU. 



s 


S 


U 


s 


U 




01, 02 


01 


01 


12 


02 


12 


12, 23 


02 


01 


12 


02 


12 


11, 22 


11 


11 


22 


11 


22 


23 


12 


02 


12 


02 


12 



During the last two time steps, each petal contributes exactly 2 particles to the central site, 
and the central site topples once. Thus after two time steps not only the petals, but also 
the central site is in the same state. Therefore, the period becomes 2. 

Next we consider SSUSS. At this stage each petal is in state 01 or 11, so if there were yet 
another S, the sandpile would be periodic with period 1. So we see SSUSSU, and because 
SSUU is forbidden, we conclude that we see SSUSSUS. 



s 


S 


U 


S 


S 


U 


s 


01, 02 


01 


01 


12 


02 


01 


12 


12, 23 


02 


01 


12 


02 


01 


12 


11, 22 


11 


11 


22 


11 


11 


22 


23 


12 


02 


12 


02 


01 


12 



At the time of the third S, each petal is in state 12 or 22. Between the third S and the fifth 
S, each petal contributes exactly two particles to the central site and returns to the same 
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state, while the central site topples once. Thus the configuration is periodic with period 
3. □ 

We conclude from the above case analysis that the activity pa is always one of 0, 1/3, 
1/2, 2/3, or 1. Table m summarizes the behavior of the periodic sandpile states for different 
values of pa- 





periodic sandpile states 


activity pa 





1/3 


1/2 


2/3 


1 


central site 


S 


s s u 


S U 


S U 


U 


U 


etals 


01 


12 02 01 


12 02 


23 12 


13 


> 22 


11 


22 11 11 


22 11 


22 11 


22 






00 




13 12 









Table 4. Behavior of the central site and petals as a function of the activity p, 



The following theorem shows that parallel chip-firing on the fiower graph exhibits four 
distinct phase transitions where the activity pa jumps in value: For each a e {0, |, |, |, 1}, 
there is a nonvanishing interval of initial densities A where pa = a asymptotically almost 
surely. 

Theorem 5.4. Let be the unique root c»/ § + = and let Cc be the unique positive 

root off - |e"3f = (. (Numerically, C = 1-6688976 ... and C = 3.3333182 . . . .} With 
probability tending to 1 as n oo, the activity density pa of parallel chip- firing on the flower 
graph F„ is given by 



'0, 


if < A < Cc 


1/3, 


if Cc < A < 2 


< 1/2, 


if 2 < A < 3 


2/3, 


if 3 < A < C 


.1 


ifC<A. 



Proof. In a given petal, let X denote the difference modulo 3 of the number of particles on 
the two sites of the petal. Observe that X is unaffected by toppling. Let Z denote the 
number of petals for which X = 0, and R denote the total number of particles, in a given 
initial configuration. Using Table HJ we can relate Z, i?, and the activity pa- 

When Pa = 0, we have less than 2n particles at the central site, at most two particles 
for the Z petals of type X = 0, and exactly one particle for the other n — Z petals, so 
< R<2n + 2Z + {n- Z) = 3n + Z. 

When Pa = 1/3, hj considering the U time step, we have R > 2n + 2Z + {n — Z) = 3n-\-Z. 
By considering the preceding S time step, we have R < 2n + 2Z + 2{n — Z) = An. 

When Pa = 1/2, by considering the U time step, we have R > 4n, and by considering the 
S time step, we get R < 2n + 4Z + 4(n — Z) = 6n. 

When Pa = 2/3, by considering the second U step, we have R > 2n + AZ + 4{n — Z) = Qn. 
By considering the S time step, we have R < 2n + 4Z + 5{n — Z) = 7n — Z . 
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When Pa = 1, we have R > 2n + iZ + 5{n - Z) = 7n - Z. 

Since for given n and Z, these intervals on the values of R are disjoint, we see that the 
converse statements hold as well: the values of R and Z determine the activity pa- We 
summarize these bounds: 

fO if andonlyif 0<P<3n + Z 
1 /3 if and only if 3n + Z < R < An 
Pa = '{1/2 if and only if 4n < i? < 6n 

2/3 if and only if 6n < R < 7n — Z 
1 if and only if 7n — Z < R. 

Everything so far holds deterministically; next we use probability to estimate R and Z. 
By the weak law of large numbers, R/n — ?• 2A and Z/n ^ Pr(X = 0) in probability. Thus, 
to complete the proof it suffices to show 

Pr(X = 0) = -(1 + 2e"2^). (4) 



We can think of building the initial configuration ax by starting with the empty config- 
uration and adding particles in continuous time. Then the value of X for a single petal 
as particles are added is a continuous time Markov chain on the state space {0,±1} with 
transitions — )■ ±1 at rate 2, and ±1 — )■ and ±1 — )■ ±1 each at rate 1. Starting in state 0, 
after running this chain for time A we obtain 



Pt{X = 0) 
Pr(X ^ 0) 



exp{2A(P-/)} 



1 




where P 



1/2 

1 1/2 



and / is the 2x2 identity matrix. The eigenvalues of P — I are 

2 J aina V2 — [_1J- uiiiv^c 



and — |, with corresponding eigenvectors fi = [2I and V2 = \ Since m — 1 ^ 
obtain (H. 



3^1 + 3f2, we 

□ 



The following theorem describes a phase transition in the driven sandpile dynamics on the 
flower graph analogous to [Theorem 3. II for the bracelet graph. We remark on one interesting 
difference between the two transitions: for A > Co the final density p(A) is increasing in A 
for the bracelet, and decreasing in A for the fiower graph. 

For A > 0, let a\ be the configuration with Poisson(A) particles independently on each 
site of Fn- Let r]x = {cr\)° be the stabilization of a\, and let 

n-l 



a;=l 



be the final density. 



Theorem 5.5. For the flower graph with n petals, in the limit n ^ 00 we have 

(1) The threshold density (c is the unique positive root o/^ = | + \e~^^ 

(2) The stationary density (s is 5/3. 
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Figure 5. Density p(A) of the final stable configuration as a function of initial 
density A on the flower graph F„ for large n. A second-order phase transition 
occurs at A = Cc- Beyond this transition, the density of the driven sandpile 
decreases with A. 




(3) Pn(A) — p(A) in 'probability, where 

p(A) = min (^A, ^ + ^e"^ 

Proof. Part 1 follows from [Theorem 5.^ 

For Part 2, we use the burning algorithm. In all recurrent configurations on F„, the central 
site has either 2n — 1 or 2n — 2 particles. All other sites have at most one particle, and in 
each petal (except the sink) there is at least one particle. For each petal that is not the sink, 
there are two possible configurations with 1 particle, and one with 2 particles. Each of these 
occurs with equal probability in the stationary state, so the expected number of particles in 
the petals is (n — 1)(| ■ 1 + | • 2) = ^ + 0(1) as n — )• oo. Therefore, the total density is 

C = lim„^o,H^t5^ + o(l)=5/3. 

For part 3, for the driven dissipative sandpile on Fn, we first stabilize all the petals, 
then topple the center site if it is unstable, then stabilize all the petals, and so on. For 
each toppling of the center site, the sandpile loses 0(1) particles to the sink. If the center 
topples at least once, then each petal will be in one of the states 11, 01, or 10, after which 
the number of particles at the center site is R — n — Z + 0{1). Recall from the proof of 
Theorem 15.41 that R/n — )■ 2A and Z/n — )■ ^(1 + 2e~^^) in probability. Thus if A < (c, then 
~^ ~ 1 + 1^ '^'^ < 1 in probability, so the sandpile does not lose a macroscopic 
amount of sand, and p„(A) — )■ A in probability. 

If A > (c, then the number of particles that remain after stabilization is 2n + n + Z + 0{1) . 
In this case, we have Pri(A) = + o(l) — ^ | + iii probability. □ 
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6. Sandpiles on THE Cayley tree 

Dhar and Majumdar |DM90j studied the abelian sandpile model on the Cayley tree (also 
called the Bethe lattice) with branching factor g, which has degree q + 1. Implicit in their 
formulation is that they used wired boundary conditions, i.e., where all the vertices of the 
tree at a certain large distance from a central vertex are glued together and become the sink. 
(The other common boundary condition is free boundary conditions, where all the vertices 
at a certain distance from the central vertex become leaves, and one of them becomes the 
sink. The issue of boundary conditions becomes important for trees, because in any finite 
subgraph, a constant fraction of vertices are on the boundary. This is in contrast to Z^, where 
free and wired boundary conditions lead to the same infinite volume limit. See |LP10] .) 

The finite regular wired tree Tq^„ is the ball of radius n in the infinite {q + l)-regular 
tree, with all leaves collapsed to a single vertex s. In Tg^n the vertex s serves as the sink. 
Maes, Redig and Saada |MRS02] show that the stationary measure on recurrent sandpiles 
on Tg^n has an infinite- volume limit, which is a measure on sandpiles on the infinite tree. 
Denoting this measure by Pr^, if h denotes the number of particles at a single site far from 
the boundary, then we have |DM90] 

From this formula we see that the stationary density is 

For 3-regular, 4-regular, and 5-regular trees, these values are summarized below: 




Figure 6. The Cayley trees (Bethe lattices) of degree d = 3,4,5. 
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tree 


Eg[h] 


distribution of height h of sand 


q degree 


PrJ/i = OJ Prg[/i = lJ PrJ/i = 2J PrJ/i = 3J Pvq[h = A\ 


2 3 

3 4 

4 5 


3/2 

2 
5/2 


1/12 4/12 7/12 
2/27 2/9 1/3 10/27 
81/1280 27/160 153/640 21/80 341/1280 



Large-scale simulations on Tq^n are rather impractical because the vast majority of vertices 
are near the boundary. Consequently, each simulation run produces only a small amount of 
usable data from vertices near the center. 

To experimentally measure (c for the Cayley trees, we generated large random regular 
graphs Gg^n, and used these as finite approximations of the infinite Cayley tree. We used the 
following procedure to generate random connected bipartite multigraphs of degree g + 1 on 
n vertices (n even). Let Mq be the set of edges (i, z + 1) for i = 1, 3, 5, . . . , n — 1. Then take 
the union of Mq with q additional i.i.d. perfect matchings Mi, . . . , Mq between odd and even 
vertices. Each Mj is chosen uniformly among all odd-even perfect matchings whose union 
with Mo is an n-cycle. 

Most vertices of (?„ will not be contained in any cycle smaller than \ogqn + 0{l) (see e.g., 
[Bol86j ). so these graphs are locally tree-like. For this reason, we believe that as n — t- oo the 
threshold density Ac ((?„,) will be concentrated at the threshold density of the infinite tree. 

Since the choice of multigraph affects the estimate of (c, we generated a new independent 
random multigraph for each trial. The results for random regular graphs of degree 3, 4 and 
5 are summarized in Tables El and [71 We find that for the 5-regular tree, the threshold 
density is about 2.511 rather than 2.5, for the 4- regular tree the threshold density is very 
close to but decidedly larger than 2, while for the 3- regular tree the threshold density is 
extremely close to 1.5, with a discrepancy that we were unable to measure. However, for the 
3-regular tree there is a measurable discrepancy (about 2 x 10"^) in the probability that a 
site has no particles. 
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n 


^samples 




distribution of height h of sand 


^topplings 


Pr[/i = 0] 


Pr[/i = 11 

L J 


Pr[/i = 21 

L J 


-=-n log"*^/^ n 


1048576 


2097152 


1.5004315 


0.0833326 


0.332903 


0.583764 


1.263145 


2097152 


1048576 


1.5003054 


0.0833321 


0.333031 


0.583637 


1.258046 


4194304 


524288 


1.5002161 


0.0833314 


0.333121 


0.583548 


1.253092 


8388608 


262144 


1.5001528 


0.0833311 


0.333185 


0.583484 


1.247642 


16777216 


131072 


1.5001081 


0.0833311 


0.333230 


0.583439 


1.242359 


33554432 


65536 


1.5000765 


0.0833307 


0.333262 


0.583407 


1.237317 


67108864 


32768 


1.5000540 


0.0833307 


0.333285 


0.583385 


1.232398 


134217728 


16384 


1.5000382 


0.0833307 


0.333300 


0.583369 


1.227548 


268435456 


8192 


1.5000269 


0.0833308 


0.333311 


0.583358 


1.222371 


536870912 


4096 


1.5000191 


0.0833308 


0.333319 


0.583350 


1.214431 


1073741824 


2048 


1.5000136 


0.0833307 


0.333325 


0.583344 


1.212751 


oo (stationary) 


1.5 


0.0833333 


0.333333 


0.583333 





Table 5. Data for the fixed-energy sandpile on a pseudorandom 3-regular 
graph on n nodes. Each estimate of K[h] has standard deviation less than 
7 • 10~®, and each estimate of the marginals Pr[/i = i] has standard deviation 
less than 3 • 10"^. The data for E,[h] appears to fit 3/2 -|- const/y^ very well, 
and extrapolating to n — )■ oo it appears that E[/i] — >■ 1.500000 to six decimal 
places. However, apparently Fr[h = 0] — )■ 0.083331 < 1/12. 



n 


#samples 


E[h] 


distribution of height h of sand 


#topplings 


Pr[h = 0] 


Pr[/i = 1] 


Pr[/i = 2] 


Pr[/i = 3] 


-^n log^/^ n 


1048576 


2097152 


2.001109 


0.073884 


0.221887 


0.333466 


0.370763 


0.623322 


2097152 


1048576 


2.000853 


0.073881 


0.221978 


0.333547 


0.370593 


0.618848 


4194304 


524288 


2.000688 


0.073880 


0.222037 


0.333599 


0.370484 


0.620894 


8388608 


262144 


2.000584 


0.073878 


0.222075 


0.333631 


0.370416 


0.631324 


16777216 


131072 


2.000518 


0.073877 


0.222100 


0.333651 


0.370372 


0.649328 


33554432 


65536 


2.000477 


0.073877 


0.222114 


0.333664 


0.370345 


0.670838 


67108864 


32768 


2.000451 


0.073877 


0.222123 


0.333673 


0.370328 


0.691040 


134217728 


16384 


2.000434 


0.073876 


0.222130 


0.333678 


0.370316 


0.699706 


268435456 


8192 


2.000424 


0.073876 


0.222134 


0.333681 


0.370310 


0.695065 


536870912 


4096 


2.000417 


0.073876 


0.222136 


0.333683 


0.370305 


0.684507 


1073741824 


2048 


2.000413 


0.073876 


0.222138 


0.333684 


0.370303 


0.673061 


oo (stationary) 


2 


0.074074 


0.222222 


0.333333 


0.370370 





Table 6. Data for the fixed-energy sandpile on a pseudorandom 4- regular 
graph on n nodes. Each estimate of E[/i] and of the marginals Pr[/i = i] has 
standard deviation less than 3 • 10~^. 
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n 


^samples 


E[h] 


distribution of height h of sand 


#topphngs 


Pr[/i = OJ 


Pr[h = IJ 


Pr[h = 2\ 


PT[h = 3J 


Pi[h = 4J 


-irn 


1048576 


1048576 


2.512106 


0.062271 


0.166547 


0.237230 


0.264711 


0.269242 


1.666086 


2097152 


524288 


2.511947 


0.062269 


0.166579 


0.237256 


0.264727 


0.269169 


1.666244 


4194304 


262144 


2.511847 


0.062268 


0.166599 


0.237272 


0.264737 


0.269123 


1.666404 


8388608 


131072 


2.511781 


0.062267 


0.166613 


0.237283 


0.264743 


0.269093 


1.666589 


16777216 


65536 


2.511743 


0.062267 


0.166621 


0.237289 


0.264748 


0.269075 


1.667322 


33554432 


65536 


2.511716 


0.062267 


0.166627 


0.237293 


0.264750 


0.269063 


1.668196 


67108864 


32768 


2.511700 


0.062267 


0.166630 


0.237296 


0.264752 


0.269056 


1.669392 


134217728 


16384 


2.511689 


0.062267 


0.166632 


0.237297 


0.264755 


0.269050 


1.671613 


268435456 


8192 


2.511683 


0.062266 


0.166634 


0.237299 


0.264753 


0.269048 


1.675479 


536870912 


4096 


2.511680 


0.062267 


0.166633 


0.237300 


0.264755 


0.269045 


1.677092 


1073741824 


2048 


2.511677 


0.062266 


0.166634 


0.237300 


0.264755 


0.269044 


1.688093 


oo (stationary) 


2.5 


0.063281 


0.168750 


0.239063 


0.262500 


0.266406 





Table 7. Data for the fixed-energy sandpile 
graph on n nodes. Each estimate of E[/i] and 
standard deviation less than 2 • 10~^. 



on a pseudorandom 5-regular 
of the marginals Pr[h = i] has 
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7. SaNDPILES on THE LADDER GRAPH 

The examples in previous sections suggest that the density conjecture can fail for (at least) 
two distinct reasons: local toppling invariants, and boundary effects. A toppling invariant 
for a graph G is a function / defined on sandpile configurations on G which is unchanged 
by performing topplings; that is 

f{a) = f{a + A,) 

for any sandpile a and any column vector of the Laplacian of G. Examples we have seen 
are 

/(cr) = a{x) mod 2 

where x is any vertex of the bracelet graph and 

/(cr) = ^(.Ti) — a{x2) mod 3 

where Xi,X2 are the two vertices comprising any petal on the flower graph F„. Both of these 
toppling invariants are local in the sense that they depend only on a bounded number of 
vertices as n — )■ oo. 

The Cayley tree has no local toppling invariants, but the large number of sinks, comparable 
to the total number of vertices, produce a large boundary effect. The density conjecture fails 
even more dramatically on the complete graph (ITheorem 4.11) . One might guess that this is 
due to the high degree of interconnectedness, which causes boundary effects from the sink to 
persist as n — )■ oo. A good candidate for a graph G satisfying the density conjecture, then, 
should have 

• no local toppling invariants, 

• most vertices far from the sink. 

The best candidate graphs G should be essentially one-dimensional, so that the sink is well 
insulated from the bulk of the graph, keeping boundary effects to a minimum. Indeed, the 
only graph known to satisfy the density conjecture is the infinite path Z. 

Jarai and Lyons |JL07j study sandpiles on graphs of the form G x P„, where G is a finite 
connected graph and P„ is the path of length n, with the endpoints serving as sinks. The 
simplest such graphs that are not paths are obtained when G = Pi has two vertices and 
one edge. These graphs are a good candidate for = Cs, for the reasons described above. 
Nevertheless, we find that while Cc and (s are very close, they appear to be different. 

First we calculate Cs- Jarai and Lyons |JL07[ section 5] define recurrent configurations as 
Markov chains on the state space 

X = |(3, 3), (3, 2), (2, 3), (3, 1), (1, 3), (3:21, (2:3)} 



Figure 7. The ladder graph. 
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describing the possible transitions from one rung of the ladder to the next. States (i, j) and 
both represent rungs whose left vertex has i — 1 particles and whose right vertex has 
j — 1 particles. The distinction between states (3, 2) and (3, 2) lies only in which transitions 
are allowed. The adjacency matrix describing the allowable transitions is given by 

/ 1 1 1 1 1 \ 



A 



1 
1 
1 
1 
1 
V 1 




















1 



1 / 



UiVi/Z, where 



Its largest eigenvalue is 2 + \/3, and the corresponding left and right eigenvectors are 

M = (1 + 73, 1 + 73, 1 + 73, 1, 1, 1, 1) 

w = (3 + 73, 1 + 73, 1 + 73, 1 + 73, 1 + 73, 1, 1)^ 

By the Parry formula |Par64] . the stationary probabilities are given by p{i) 
Z is a normalizing constant. So 

p(3,3) 
p(2,3)=p(3,2) 
p(l,3)=p(3,l) 
p(2;3)=p(3:2) 
where 

Z= (l + 73)(3 + 73) + 2(1 + a/3)^ + 2(1 + 73) + 2. 
Thus we find that for the ladder graph in stationarity, the number h of particles at a site 
satisfies 



(l + 73)(3 + 73)/Z 

(1 + 73)Vz 

(1 + 73) /z 
1/Z 



PT[h 


= 0] = 




1 1 

2 ~'~ 3 


= 0.0773503 


Fr[h 


= 1] = 


5 
4 


7^3 _ 
12 


= 0.2396370 . 


Fr[h 


= 2] = 


1 
4 


_A_ Vs _ 
4 


0.6830127.. 




E[h] = 


7 
4 


_ 

12 


1.60566243 



In contrast, the threshold density for ladders appears to be about 1.6082. Table |8] summarizes 
simulation data on finite 2 x n ladders. 
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n 


7f samples 


K[h\ 


distribution of height h of sand 


#topphngs 


Pr[/i = 0] 


Pv[h = 1] 


Pr[/i = 2] 




256 


4194304 


1.60567 


0.07695 


0.24043 


0.68262 


0.094773 


512 


2097152 


1.60693 


0.07656 


0.23996 


0.68349 


0.095366 


1024 


1048576 


1.60757 


0.07636 


0.23970 


0.68393 


0.095864 


2048 


524288 


1.60788 


0.07626 


0.23960 


0.68414 


0.096316 


4096 


262144 


1.60805 


0.07621 


0.23952 


0.68426 


0.096545 


8192 


131072 


1.60814 


0.07618 


0.23950 


0.68432 


0.096753 


16384 


65536 


1.60816 


0.07618 


0.23949 


0.68434 


0.097113 


32768 


32768 


1.60818 


0.07617 


0.23948 


0.68435 


0.096944 


65536 


16384 


1.60820 


0.07616 


0.23948 


0.68436 


0.097342 


131072 


8192 


1.60820 


0.07617 


0.23946 


0.68437 


0.097648 


262144 


4096 


1.60821 


0.07615 


0.23949 


0.68436 


0.096158 


oo (stationary) 


1.60566 


0.07735 


0.23964 


0.68301 





Table 8. Data for the fixed-energy sandpile on 2 x ladder graphs. Each 
estimate of E,[h] and of the marginals Pr[/i = i] has a standard deviation 
smaller than 10~^. To four decimal places, the threshold density (c equals 
1.6082, which exceeds the stationary density (s = 7/4 — \/3/12 = 1.6057. The 
total number of topplings appears to scale as n^^^. 
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